R-matrix propagation with adiabatic bases for the photoionization spectra of atoms in 

magnetic fields 

F. Mota-FurtadcQ and P. F. O'MahonjQ 

Department of Mathematics, Royal Holloway, University of London, Egham, Surrey TW20 OEX, United Kingdom 

(Dated: February 1, 2008) 

The photoionization spectrum of an atom in a magnetic field is calculated by combining R-matrix 
propagation with local adiabatic basis expansions. This approach considerably increases the speed 
and the energy range over which calculations can be performed compared to previous methods, 
allowing one to obtain accurate partial and total cross sections over an extended energy range for 

■ an arbitrary magnetic field strength. In addition, the cross sections for all atoms of interest can be 
calculated simultaneously in a single calculation. Multichannel quantum defect theory allows for a 

, detailed analysis of the resonance structure in the continuum. Calculated cross sections for a range 

■ of atoms in both laboratory and astrophysical field strengths are presented. 
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f^h I. INTRODUCTION 

Ph. 



The spectrum of an atom in a magnetic field has played an important role in the development of quantum theory 
and in atomic structure For a long time the effect of the external field on the atom was treated perturbatively. 
However with the discovery of large magnetic fields in white dwarf stars ( 10 2 — 10 5 T) and neutron stars ( 10 7 — 10 9 T) 
in the seventies attention has focussed on non-perturbativc treatments of the field-atom interaction . 2']. For laboratory 
strength magnetic fields ( 6T) the potential due to the applied field only becomes comparable to the intrinsic Coulomb 
potential for an electron in a high Rydberg state or continuum state of the atom. In the eighties and nineties research 
was focussed on atoms in laboratory strength fields as these systems have an inherent non separability arising from 
the competing spherical symmetry of the atom and the cylindrical symmetry of the applied field which leads to 
O ' the classical system exhibiting chaotic behavior Q. An atom in a field thus provided an experimentally realizable 
1—1 1 quantum system whose corresponding classical phase space is chaotic, serving as a prototype for studying classical 
and quantum chaos. This lead to fruitful developments in the theory of quantum chaos [4(. 

For bound states of atoms in a magnetic field large scale basis set calculations have proved very successful in finding 
ly-) ' the energy eigenstates and photoabsorbtion spectra of atoms in moderately strong fields 0. For super strong fields 
, such as those found in neutron stars, the field can modify the atomic structure of the ground state of the atom and 
different theoretical techniques have to be used. Currently for these cases the energy eigenstates and photoabsorbtion 
spectra are only known for some low levels of a few light atoms j5|]. 
Q\ • The positive energy or continuum spectra of an atom in a field proved more challenging particularly in calculating 
photoionization cross sections at laboratory strength fields. The three main theoretical methods that have been 
successful at calculating cross sections at both laboratory and astrophysical strength fields are the complex coordinate 
method of Delande et al @, the R matrix method of O'Mahony and Mota-Furtado @, Q, and the diabatic by sector 
method of Watanabe and Komine Q . A detailed comparison between theory and experiment has been possible due 
to the high resolution experiments carried out by Iu et al 10] on lithium in a field of about 6T. Although all three 
approaches have recreated the experimental spectrum of Iu et al over a narrow energy region they are not particularly 
suited to calculating the photoionization spectrum over a large energy region. We present here a major improvement 
on previous R matrix methods applied to this problem by using local adiabatic basis states to propagate the R-matrix 
from low r to the asymptotic region. This leads to a large saving in both the CPU time and computer memory 
required to perform the calculation. In addition, for a given value of the field strength, the cross sections for all atoms 
of interest can be calculated in one step without the need for any additional propagations. We demonstrate that 
by using a combination of R matrix propagation with local adiabatic basis states and multichannel quantum defect 
theory (MQDT) [ll[ one has a method that can be used to calculate photoionization cross sections of any atom over 
very large energy regions and field strengths. MQDT can also be used to analyze the resonance structure in detail. 
An efficient approach to calculating such cross sections is of importance in many areas where magnetic fields play a 
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role, for example in calculating stellar opacities for magnetic white dwarfs or for evaluating recombination rates for an 
atom or ion in a magnetic field at low temperatures (e.g. anti-hydrogen) where one has to calculate the cross section 
over very large energy regions for a given field strength. 

In section|TI]we give the theory used to evaluate the photoionization cross section of an atom in an external magnetic 
field by combining R matrix propagation with local adiabatic basis states and MQDT. In section IIIII all the required 
details of the computation are given. The results of the calculations are presented in section HVl for a variety of atoms 
at both laboratory and astrophysical field strengths and a concluding section is given in section |Vj In the Appendix 
we give details on how to construct the Hamiltonian matrix for the propagation. 

II. THEORY 

The Hamiltonian for a hydrogen atom in a magnetic field (taken to be in the z direction) in the symmetric gauge 
can be written using atomic units (h = m = e = 1) as 

H = -^-i+/3L 2 + i/3Vsin 2 f? (1) 

where the magnetic field B is measured in atomic units by (3 = B/Bq with Bq = 4.70108 x 10 5 T. (We neglect the spin 
as it only produces a uniform shift in the energy scale) . The Hamiltonian has two conserved quantities in addition to 
the energy, namely the z-component of the angular momentum L z and tt z the z-parity. The eigensolutions can thus 
be studied for fixed values of m the azimuthal quantum number and for tt z = ±1. The linear Zeeman term in eqn. 
(TT|) thus only adds a uniform shift to the total energy. 

When a photon excites an electron to the continuum, for a given magnetic field strength, one can in general identify 
three regions of interaction with the continuum electron |8[. Typically for low r the spherically symmetric Coulomb 
potential dominates over the cylindrically symmetric diamagnetic term or quadratic term in eqn. |T]), at intermediate 
values of r the Coulomb and magnetic potentials are of comparable strength (the strong mixing region) and at high 
values of r or asymptotically in r the cylindrical symmetry of the diamagnetic potential predominates. For a general 
atom (or molecule) one adds a fourth region, the core, where the excited electron interacts with the multi-electron 
core before emerging into the Coulomb region described by the Hamiltonian in eqn. ([1]). Exploiting this natural 
partition in configuration space forms the basis of the R-matrix approach to solving atomic and molecular problems 
where solutions are sought in each region and then matched together at the boundaries between the regions to form 
the solution over all space [12]. A novel aspect of the magnetic field problem is having to deal with the change 
in symmetry from spherical to cylindrical which involves introducing two-dimensional matching procedures [§|. We 
describe below how the R-matrix is propagated through the regions described above and how the two dimensional 
matching procedure is implemented to give the reactance matrix and the photoionization cross section. 

A. Propagating the R-matrix 

An atom in a magnetic field is assumed to be excited from an initial state, either a ground or low lying excited 
state, by a polarized photon leading to an electron in the continuum with specific values of m and n z . The electron 
emerges from the first region, the core region, into the second or Coulomb region with energy e where eqn. ([1]) can be 
approximated by the field free Hamiltonian because the diamagnetic terms are negligible in comparison with those 
from the Coulomb potential of the atomic core. (We shall assume here that the field strengths for non-hydrogenic 
atoms are not large enough to significantly distort the core (i.e. (3 < 1). In this case a different treatment would be 
required for the core although the propagation outside the core could still be implemented) . Therefore at some radius 
r = a in the Coulomb region the radial form of the wavefunction, Ff(r), of the continuum electron can be written in 
terms of a linear combination of the energy normalized regular s and irregular c Coulomb functions [l3| in spherical 
coordinates or a phase shifted Coulomb function giving the general form of the solution as 

* e = mr)Yim(0, <f>)=J2 A i ( s t( r ) + c i W tan (^)) Y lm {6, cf) (2) 
i i ^ ' 

where the are constants to be determined. The quantum defects /z; represent the effects of the non-hydrogenic field 
free core 11] and can be calculated ab-initio or taken from experiment. (The more general case of a multi-channel 
wavefunction with a reactance matrix instead of quantum defects in eqn. ^ is straightforward to include). Yi m (0, <p) 
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are the spherical harmonics. Knowing the phase shifted Coulomb function and its derivative at r = a, the R matrix 
[l2| can be constructed on the outer boundary of this region as 



-i -l 

Sw, (3) 

since the coefficients A\ cancel out in the above expression. 

In the third or strong mixing region the effects due to the Coulomb potential of the core and those due to the 
magnetic field are of a comparable size. This region is defined by a < r < b where the radius b is taken to be large 
enough such that the Hamiltonian is separable in cylindrical coordinates. The change in symmetry of the potential, 
from spherical to cylindrical, is therefore completely contained within this region. We wish to propagate the initial 
R-matrix in eqn. ([3]) to obtain the R matrix at the outer boundary of this region r = b [lil Il5| |. To do this we 
divide the region into N radial sectors with radii a — > a\, a\ — > a^, • • • , ojv-i — > b. The size of each sector and the 
number of sectors N are important parameters in the calculation and we show how these are optimally determined 
later. Within each of the sectors we construct a local adiabatic basis as follows. For the n th sector we take a radius 
r™ within the sector a n -\ < r" < a n , which is usually the mid point of the sector, and we diagonalize the fixed r or 
adiabatic Hamiltonian H a d 

H ad K; 9, 0) = -i^ - 1 + \p 2 Kf sin 2 6 (4) 
2(C) C 2 

in a basis set of spherical harmonics such that 

&(r?;M)=X;dlAljm(M), (5) 
i 

where the di\ are constants giving the eigenvalue equation 

H ad (j>x = U x (r2)4>x- (6) 

The adiabatic states form a locally optimized basis set for each sector and the adiabatic potential curves can be 
produced by plotting the eigenvalues obtained from the diagonalization at successive r™. A typical set of curves is 
shown in figure [T] At small r the adiabatic eigenfunctions functions are spherical harmonics and the potential exhibits 
the centrifugal barrier. At large r these curves have the equal energy spacing of Landau states (combined with the 
— 1/r fall off) and the eigenfunctions arc localized in the cylindrical coordinates p and <f>. The change in symmetry 
happens predominantly around a region of avoided crossings that can clearly be seen in the diagram. The functions 
4>\ are therefore a very good basis with which to represent the angular part of the wavefunction in the local region 
around r™, namely within the sector n. 

To propagate the R-matrix from sector to sector we use these basis states within each sector and we diagonalize 
the full Hamiltonian (eqn. Jl}) plus the Bloch operator L or surface term [l6j . 



Riv 



s\ (a) + c\{a) tan(7r/i;) 



\ sK r )+c e i( r ) tan(TT^) 



L = ^(s(r-a n )^-5(r-a n -i)-^j, (7) 

in a basis set consisting of a product of orthogonal radial functions fj(r) and the adiabatic functions generated for 
that sector (f>\(r2;9,4>)- The radial basis functions used are defined in terms of Legendre polynomials Pj as follows 

El 



m = \ 2j 1 Pi-iiv), (8) 

V a n — o,n-l 



where u = 



a n + a n -i 



071— 1 y y 2 
The eigenvalue equation is thus 



(H + L)^ k = E k ^ k . 



(9) 



Radius r , (a.u.) 



6000 



FIG. 1: The first 20 adiabatic eigenvalue curves obtained from diagonalizing the adiabatic Hamiltonian H a d in a basis set of 
spherical harmonics at successive r. The magnetic field strength used was 6T. The centrifugal barrier can be seem at small r 
and the Coulomb plus equal Landau spacing is seen as r — > oo. There is a set of avoided crossings in between. 



The eigenfunctions obtained from diagonalizing this operator are therefore 

* fc = £4^0 A (C;0,# (10) 

The total continuum wavefunction at any energy e, , 5 e , can be expanded in terms of these R-matrix eigenstates 1 $> k . 
Since the general solution in the n th sector can also be written as \l/ e = J2\ P\{ r ) ( t>\{ r a^ ®i 4 1 ) it is straightforward to 
show [l5j], using the operator H + L and eqn. ([9]), that the values of the functions and derivatives on the boundaries 
of the sector are related by 



F(on-i) = r 2 F'{a n ) - nF'(a„_i) 
F{a n ) = r 4 F'(a n ) - r 3 F'K-i), (11) 

where the matrix elements of ri to r^, called the sector R- matrices, are given by 

t n\ _ 1 gik(a n -i) gjk(a n -i) _ 1 9ik(a n -i) 9jk(a n ) 

[lHj 2^ E k -e [T2h3 2^ Ek-e 



/ n\ _ 1 9ik{a n ) 9jk(a n -i) ,4^ _ 1 9ik{a n ) 9jkWn) 
\3Jio 2 2^ E k -e [ iHj 2^ E h -<= ' 

k k 

and 



9\k 



« = E^^- (13) 



.1 



In summary, knowing the eigenvalues and eigenvectors of eqn. ([9]) one can construct the sector R-matrices ri to 
r4 above which, through eqns. (fTTj) . relate the radial solutions and their derivatives on the boundaries of the sector. 
The R-matrix relates the function to its derivative, i.e. F(a n ) — R(a n )F'(a n ), and a simple manipulation of eqn. 
(fTTj) yields the relationship between the R-matrix on the inner and outer boundaries of the sector 



R(a n ) = r? - r 3 " I r? + R(a„_i) 1 r£ (14) 
where R(a„) and R(a n _i) are represented in the same adiabatic basis set as the sector R matrices. 
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As the adiabatic basis changes from sector to sector we need finally to change the basis representation of the R 
matrix. The matrix with elements 

{T n -^) xy = (0 A (rrV,0)l0v(C;M)> (15) 
is thus constructed. One uses this transformation to change the basis representation of R giving 

R = T B_1 ' n Br- 1 ' B ) (16) 



where T is the transpose of T. 

Starting with some initial R matrix it can thus be propagated from sector to sector using eqns. (fTl)) and (|16[) . 
Using the initial input R-matrix given by eqn. ([3]) the propagation gives the final R matrix on the outer boundary in 
the asymptotic region at r = b with the final R matrix being represented in the local adiabatic basis of the last sector. 

Although it is possible to propagate the R matrix itself at each of the sector radii as described above, it is more 
practical and efficient to derive global sector R matrices (Ri,R,2,R.3 and R4) relating the first and n th sectors 
which can be built up sequentially using eqn. (fTTjl . One initially generalizes eqn. (| 1 1 1) to relate the values of the 
functions and derivatives on the boundaries of the first and n th sectors 

F(a x ) = TL%F'(a n ) - R?F'(a x ) 

F(a n ) = RJF'K) - Rjj"(ai). (17) 

The operator relations for these global sector R matrices, including the change in basis, have been derived by Stechel 
et al [171] and can be obtained by matching the wavefunction and its derivative on the boundaries between each of the 
n sectors . They are 

On 1 tt> n — lriT^ -^ryn rr-m — l,nT) n — 1 

1 1 — I V — 1 \ , 1 A 1 III 

R£ = R' 2 l - 1 T Il_1 '"Z"r™ 
■pji r 7l Z n T Tl— 1 

Rn _n nrjn n 
4 — r 4 — r 3 A r 2 



where Z™ = 



^.n + r -l,» R jf ^ . ( 18 ) 



In these equations R™ are the global sector R matrices for all sectors up to n, r" are the sector R matrices for sector 
n and T n_1: " is the transformation matrix between the adiabatic basis used in sectors n — 1 and n. In an analogous 
way to eqn. (| 14[) the global sector R matrices at the end of the propagation through all N sectors are used to relate 
the R matrix at r — a to the R matrix at r = b 

R(6) = Rf -Rf (Rf + R(a))~ 1 R 2 v . (19) 

Note that R^ to R^, determined from eqn. ([T5]) above, are independent of which atom one uses so that once they 
are calculated from the propagation then R(6) for a whole set of atoms can be evaluated at once using eqn. (Til?)) 
by just using the appropriate quantum defects to calculate R(a) in eqn. ([3]). Hence for given values of B, m and 
7r 2 , the cross sections of all atoms of interest can be calculated simultaneously without the need for any additional 
propagations. 



B. Asymptotic region 

Having found the R-matrix at r — b by propagation we must match this to the asymptotic solutions at r = b to 
find the solution over all space. For large r, r > b, the magnetic field dominates. Since the motion in p is bounded 
— — * — , and the Hamiltonian in eqn. ([T]) is separable in cylindrical coordinates 

ld?_ _ 1 „ „ ( 1 
2d? ~ ~ 



H = —^--+ H L + ° 3) ( 20 ) 



where is the Hamiltonian for the Landau states 

119/ 9 \ m 2 I 
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Hl has eigenvalues Ef = (2i + \rn\ +m+ 1)0, i = 0, 1, 2 . . . and its eigenfunctions are the Landau states $>i(p, 4>) El- 
The asymptotic region c<z<oo 7 0<p<oo with c less than the radius r — b, is therefore chosen to conform with 
the cylindrical symmetry of the problem. In this region a set of j linearly independent solutions may be written, as 
is standard in scattering theory [ll| , in terms of the solutions of eqn. (j2"0|) . These are a product of Landau states <E>i 
and a linear combination of energy normalized regular and irregular Coulomb functions in z, s and c, evaluated at an 
energy a — e — Ef', namely, 



(22) 



The constants Kkj, the reactance matrix or K matrix, are to be determined by the matching procedure. 

The R matrix R(6) having being evaluated using eqn. (|19|) is now matched through a two dimensional matching 
procedure to these asymptotic solutions on an arc at r = b 8] . To perform this matching the integrals 



b x v ej dn 



(23) 



must be evaluated, i.e. the asymptotic solutions must be projected onto the local adiabatic solutions on the radial 
arc at r = b. This is done by evaluating numerically the four matrices P, Q, P', Q' with elements 



Pxj(b) = J 


Mb; 6 


i 


) Sy'O) 


dn 

r—b 


Qxj(b) = J 


MM 


i 


) Cij(z) 


dn 

r—b 


P'xj(b) = J 


Mb; 6 


i 


) s'ijiz) 


dn 

r—b 


Qxj{b) = J 


Mb;t 


i 


) %(z) 


dn 

r—b 



(24) 



where ' indicates the derivative with respect to z. This gives the regular and irregular components of the solutions at 
r = b and allows one to calculate the outer R-matrix from the asymptotic region at r = b in terms of K the reactance 
matrix. Equating the inner and outer R-matrices at r = b gives 



R 



(Q K) 



P' + (Q' K) 



Re-arranging this expression gives us the equation for the K matrix 



K 



n -1 



(R Q ) Q 



(R P ) P 



(25) 



(26) 



Knowing K we have the energy normalized solution over all space and we can calculate both partial and total 
photoionization cross sections. 



1. Multichannel quantum defect theory 

For a given total energy e there are two possibilities for the behaviour of the Coulomb solutions in eqn. (f2"2")l . If 
€i = e — Ef > the channel is open and the Coulomb functions oscillate at r = b and all the way to infinity. If 
ti < then the channel is closed and the solutions must decay as z — > oo. The physical K matrix therefore is a square 
matrix with the dimension of the number of open channels. However with multichannel quantum defect theory one 
exploits the known analytic properties of the Coulomb functions to enforce the boundary conditions for the closed 
channels. For the closed channels there are two possible scenarios at r = b. Either the channel is strongly closed 
and is already exponentially small at r = b or it is weakly closed, i.e. the Coulomb functions are still oscillating at 
r = b before decaying at infinity. In MQDT, for a weakly closed channel, one instead uses Coulomb functions in eqn. 
(I22p which don't decay at infinity and one treats the channel as if it is open only enforcing the boundary condition 
at infinity analytically in a final step [ll|. Therefore the resonance structure due to these weakly closed channels 
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can be calculated analytically. Hence the K matrix calculated by doing the matching using MQDT, denoted by /C, 
has dimension given by the number of open plus weakly closed channels. In general K. has a smooth dependence on 
energy because the Rydberg series of resonances converging on the Landau thresholds corresponding to the weakly 
closed channels have not as yet been included. 

The open part of the actual physical reactance matrix K can be recovered from the matrix K, by the formula [Til ] 

K DO = K, 00 - K, oc ^tan(7w) + JC c< ^j IC co . (27) 

The o and c subscripts refer to the open and closed channels of K, and the tan7rz/s form a diagonal matrix where the 
i/s are related to the energies a by 

Mi = - (28) 

This way of obtaining the K matrix has two major advantages. Firstly it is much quicker computationally since K, 
represents a reactance matrix with a lot of the resonance structure removed it varies much more slowly with energy 
than the full reactance matrix K. This allows one to calculate K, on a coarse energy mesh with fairly large energy 
spacings. These /C's can then be used to calculate K over an arbitrarily fine energy mesh using the analytic formula 
in eqn. (|27p . This allows the propagation stage to be performed at fewer energy points ultimately speeding up the 
calculation enormously. The second major advantage of this approach is that the resonance structure converging to 
a particular Landau threshold can be identified by removing it from K. This is done by keeping open the relevant 
Landau channel in the evaluation of equation (|2T[) . By comparing a spectrum with all resonances converging on 
a particular Landau threshold removed with that of the full spectrum it is possible to determine which resonances 
converge to which thresholds [1, [l8| . This technique is demonstrated in section HVl 



C. Photoionization cross section 

The photoionization cross section is given by 

a = 47r 2 aw|(*7|e.r1* )| 2 (29) 

where a is the fine structure constant, oj the photon ener gy, e the polarization direction, ty the initial bound state 
and the energy normalized 'incoming' wavefunction [191 ]. There is a standard transformation to go from the 
K-matrix form given in (|22p to the S-matrix or incoming form ^~ [ll| . Once the K matrix and hence the asymptotic 
form is known from the matching, the photoionization cross section is evaluated by calculating the amplitude of the 
wavefunction near the origin, and hence the dipole integrals, as follows. 

The 'incoming' wavefunction on the inner boundary of the strong mixing region (r = a) can be written as 

*e" =£)Ff(a)15m(M)- (30) 
I 

The wavefunction on the outer boundary of the strong mixing region can be written as 

K =J2 G x( b )Mb;0^), (3i) 

A 

where is the energy normalized 'incoming' asymptotic solution calculated from the matching procedure in section 
IIIBI The radial solution F~ and G~ are linked by the global sector R matrices as given in eqn. (flT)) so that 

F-(a) = R^G'-(6)-RfF'-(a). (32) 

The radial solutions on r — a can be written as F~ = SA~ from eqn. ^ where S is a diagonal matrix with elements 

S; (a) + c\ (a) tarn TTfn (33) 

and A~ are the field and energy dependent amplitudes. Substitution into equation eqn. p2[) yields the equation 

(S + Rf S')A~ = Rf G'- (b). (34) 
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The coefficients A~ can therefore be evaluated by solving the set of linear equations because G'~(b) is known once the 
K matrix is known. Using the coefficients A" one can then simply express the cross section in terms of the field free 
photoionization cross section as the coefficients A~ give the difference in amplitudes between the field dependent and 
field free amplitudes . For example for excitation from the Is state of hydrogen using linear polarized light (Am = 0) 
only the I = 1 component of the final state will be accessed. The photoionization cross section in this case is therefore 
given by 

a(e) =\A^(e)\ 2 a B=a . (35) 
where <jb=q is the field free photoionization cross section for hydrogen. 



III. COMPUTATIONAL DETAILS 



The first point to address is the choice of the radii a and b. The inner radius must be taken larger than the core 
(~ 1 a.u.) for non-hydrogenic atoms yet small enough that the diamagnetic term is still negligible compared to —1/r. 
For example we took a — 200 atomic units for laboratory strength fields. For the asymptotic radius b, extensive 
calculations show that it is surprisingly large indicating that the long range coupling due to the terms of 0(1/ z 3 ) falls 
off very slowly. Based on our experience we used the following empirical formula for an arbitrary field 

io-^ 2/3 



b = 700 I — 1 . (36) 

The number of sectors between r = a and r = b and their sizes are chosen in the following way. First one chooses 
the maximum energy to be used in the calculation, e max . The transformation matrix (T n_1 ' n )\\> between adiabatic 
functions evaluated at two radii r™^ 1 and r™ for consecutive sectors can be calculated using eqn. (|15p. If r™ -1 was 
equal to r™, (T ,1 ~ 1,,1 )aa' would be the identity matrix. As the distance | r™ -1 — r™ | increases the corresponding 
adiabatic functions become more different. The effect this has on T is that the diagonal elements get smaller and 
the off diagonal elements get bigger. The maximum sector size for the angular adiabatic basis is found by choosing a 
limit on how small any of the diagonal elements of T can become. Choosing a value of 0.5 for the smallest diagonal 
element restricts the size of each sector for a laboratory strength field to those shown in figure [2j One can see that 
for small and large r the adiabatic functions don't vary much with r being close to spherical harmonics and Landau 
states respectively leading to large sector sizes. The intermediate range of r, where the avoided crossings in the 
potential curves shown in figure [1] are present, is where the angular functions are changing rapidly and one requires 
small sectors. 

The radial basis sector widths are chosen by comparison with the local wavelength of the Coulomb functions, namely 

constant 

Radial sector size = r . (37) 

[2 (e max + i)]* 

If a radial basis of ten Legendre polynomials is taken, it is found that a constant = 6.0 is sufficient to produce an 
accurate description of the wavefunction [1 51 ] - For the maximum energy used, the sectors resulting from this criterion 
are also shown in figure [5] There are therefore two criterion for choosing the sector sizes; one from the radial basis 
and one from the angular basis. The smallest value of these two determines the actual sector sizes and these are also 
shown in figure [U 

The adiabatic functions 4>\(t™;6,4>) are obtained by diagonalizing the adiabatic Hamiltonian H a d m a basis of 
spherical harmonics. The number of spherical harmonics used must ensure the functions (pxir^O^cj)) are properly 
converged. At the matching radius r = b, the ratio between the diamagnetic term and the potential term in H a d will 
be at its greatest. It is therefore at this radius that the largest number of spherical harmonics will be required because 
it is here that the adiabatic functions will have their greatest degree of cylindrical symmetry. We choose the number 
of Vs so that the eigenvalue corresponding to the second closed channel at r = b for an energy e ma x is accurate to 
0.5%. It has been verified that this is sufficient to give convergence in the final cross section. 

When the R matrix is propagated from r = a — ► r = b, it is necessary to retain, within any one sector, all of the 
locally open channels plus a few of the locally closed ones. The number of closed channels required is directly related 
to the threshold value on the adiabatic angular functions used to determine the sector sizes. The number of closed 
channels needed, however, is constant for all of the sectors. For a given energy e, because the number of locally open 
channels changes with radius (see figure [T]), the total number of channels retained in any one sector varies. Because it 
is very difficult to include a channel half way through the propagation, the number of channels retained in a sector is 
found by the following method. Firstly the maximum number of open channels retained for a given e max in any one 
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FIG. 2: The sector radii obtained using the angular condition only, the radial condition only and both conditions as described 
in the text. The first and last sector radii are 200 atomic units and 12000 atomic units respectively. The magnetic field strength 
is 6T. 

sector is found and this sector is labeled k. In every sector up to and including k the number of channels retained 
is this maximum number of open channels plus the number of extra closed channels. After sector k the number of 
channels retained is the number of locally open channels plus the number of closed channels. The redundant channels 
in the global sector R matrices are removed for any given energy e by a simple truncation as it is required during the 
propagation stage. Essentially once the saddle point of the adiabatic curves (see figure [T|) is passed the number of 
channels retained in the R matrix propagation can be reduced. An example of the number of channels retained in a 
program run for a magnetic field strength 470T is shown in table [J showing that from about 190 a.u. on the number 
of channels propagated decreases with r. 



Outer sector radius 


no. of local open channels 


no. of channels used 


23.41 


3 


8 


43.94 


3 


8 


72.05 


4 


8 


108.04 


5 


8 


141.02 


6 


8 


189.01 


4 


6 


247.23 


3 


5 


313.78 


2 


4 


388.78 


1 


3 


472.01 


1 


3 


563.75 


1 


3 


663.93 


1 


3 


700.00 


1 


3 



TABLE I: Table of the number of open channels and the total number of channels used in each sector of the propagation for a 
magnetic field strength of 470T. Two extra closed channels are retained in each sector. 



The number of channels needed to match the R matrix to asymptotic solutions in order to obtain the JC matrix was 
checked for each spectrum to ensure convergence. Once the K, matrix was obtained, however, the number of channels 
needed for the calculation of the cross section was reduced to the number of open channels plus the number of weakly 
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FIG. 3: Photoionization cross section of lithium in arbitrary units versus the energy in wavenumbers measured relative to the 
field free ionization threshold. Excitation is from the 3s state in a field strength of 6.1143T using linearly polarized light. The 
final state has m = and -k z = — 1. The top panel shows the full spectrum and the second, third and fourth panels show 
the spectrum with all the resonances converging on the nearest, second nearest and third nearest threshold removed using 
multichannel quantum defect theory. 

closed channels. Weakly closed channels were taken to be those that were within an energy (in a.u.) of i of their 
corresponding Landau energy. 

The cross section is obtained over a coarse mesh of energies initially and then MQDT is used to calculate the cross 
section over an arbitrarily small energy range. The convergence of the cross sections was tested by varying the inner 
and outer radii , a and 6, and the other parameters in the calculations. 

IV. PHOTOIONIZATION CROSS SECTIONS FOR LABORATORY STRENGTH FIELDS 

Calculations at laboratory strength fields are the most demanding numerically as the radial distances over which 
one needs to propagate the R-matrix are large. The number of channels needed in the calculation can also become 
quite large. We focus on the photionization spectrum of lithium in a magnetic field of 6.1143T (/3 = 1.3 x 10 -5 ) 
where experimental data exist (Iu et al [lOj ) . Lithium is excited from the 3s state using linearly polarized light giving 
a final state with m = and tt z = —1 or odd z-parity. We show in figure [3] the calculated spectrum between the 
photoionization threshold (i.e. i = or first Landau level) and the second Landau level or threshold i = 1. 

The only significant quantum defect is for I = 1 , n p = 0.053. (The I = continuum state doesn't play any role 
for excitation from s states.) The photoionization cross section is given in arbitrary units but it can be put on an 
absolute scale if the field free cross section is known. The radius a is 200 and b — 12600. The adiabatic matrix 
threshold was taken to be 0.1. Together with an e m ax of 3.9 x 10~ 5 a.u. this gives the number of sectors to be 64. The 
maximum number of locally open adiabatic channels is 27 and taking 13 extra closed channels the maximum number 
of channels overall is 40, hence the maximum size matrix to be diagonalized is 400 since there are 10 radial basis 
functions per sector. The cross section is calculated over a course mesh at 400 energy points between the thresholds. 
The full spectrum is then calculated semi-analytically using MQDT and eqn. (|27[) with over 10,000 energy points. 
The spectrum obtained from this calculation is displayed in the top panel of figure [3] showing the complete resonance 
structure between the thresholds. MQDT can also be used to examine the resonance structure by keeping individual 
weakly closed channels open, i.e. by not enforcing the closed channel conditions asymptotically. For a single threshold 
this is equivalent to a Gailitis average [HI [111 over the resonances converging to that threshold. For several thresholds 
this is a type of generalized Gailitis average over the resonance structure and would be similar to convoluting the 
actual spectrum with a Gaussian of a certain width [181 ]. The differences between the first and second panel allows 
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FIG. 4: Photoionization cross section of lithium, in arbitrary units, over a very small energy range just below the first excited 
Landau threshold. This cross section demonstrates the arbitrary resolution of the theoretical technique. 



one to identify the perturbed resonances which converge to the second Landau level or threshold i = 1. The third 
panel in figure [3] has resonances converging on the next two thresholds removed (i = 1 and i = 2) and the final 
panel has resonances converging to the i — 3 threshold also removed. Some of the remaining modulations may be 
resonances converging on higher Landau levels but cannot be removed because the wavefunctions of these resonances 
are completely contained within the radius b — 12600. A lower matching radius would be required to further reduce 
the spectrum. Note that the final panel contains only about 20 resonances compared to the very large number in the 
first panel indicating that much of the complex resonance structure is due to several Rydberg series interacting with 
a finite number of short range perturbers. 

In addition to helping in the analysis of the spectrum, MQDT allows one to calculate the resonance structure to an 
arbitrarily small resolution. This is demonstrated in figure [4] where an enlargement of the spectrum in a very small 
energy region just below the first excited Landau level is shown. 

However the real power of the combination of R matrix propagation with adiabatic bases can be seen in figure 
where we calculate the photoionization cross section over an extended energy region covering over 6 Landau thresholds 
from the ionization threshold. Because the number of open channels increases with energy the size of the matrices 
to be diagonalized are slightly larger. The total number of sectors increases also but even for the highest energy, the 
number of channels didn't exceed 50 giving the largest matrices that need to be diagonalized to be of the order of 
500. The cross section was calculated on a coarse energy grid of 500 points over each of the thresholds before applying 
MQDT. The method scales in a reasonable way thus allowing one to calculate the photoionization cross section of an 
atom over a large energy range above the ionization threshold. 



V. PHOTOIONIZATION CROSS SECTIONS FOR ASTROPHYSICAL STRENGTH FIELDS 

Photoexcitation and photoionization cross sections of light elements such as hydrogen and helium are important 
in understanding the properties of white dwarf and neutron stars 0] • The method detailed in section [TT] can equally 
well be applied to atoms in astrophysical strength magnetic fields. In fact the computation is much easier in this 
case compared to laboratory strength fields as the radius b is much smaller and the number of channels and sectors 
is smaller too. We give just two examples to illustrate the point. We first consider the photoionization spectrum of 
hydrogen from the ground state using linearly polarized light in a field of 23,500T (/3 = 0.05 a.u.) (See figure 0) The 
radius a was 1 a.u. and b — 50, the number of sectors used was 20 and up to ten channels were used in the propagation. 
This spectrum has also been calculated using the complex co-ordinate method @, [2l[ . Excellent agreement is found 
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FIG. 5: The extension of the photoionization spectrum of lithium in a magnetic field of 6.1143T shown in figure[3]to an energy 
range covering over 6 Landau thresholds. Each panel shows the spectrum over one Landau threshold. 
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FIG. 6: Photoionization cross section of hydrogen in megabarns versus the energy in a.u. measured relative to the field free 
ionization threshold. Excitation is from the Is state in a field strength of 23, 500T using linearly polarized light. The final state 
is m = and n z — — 1. The energy range covers the first couple of Landau thresholds. 



between the two methods away from the energy region near to the ionization thresholds. As the complex co-ordinate 
method uses a finite basis it cannot represent all of the Rydberg structure just below the ionization thresholds. 

The second example is the photoionization spectrum of helium from the ground state in a field of 4700T (j3 — 0.0001 
a.u.). The only significant quantum defect is for I = 1, fi p = —0.012. (The I — continuum state will only play a 
role for excitation from ls2p state.) The radius a was 5 a.u. and b = 150, the number of sectors used was 20 and 
a maximum of twelve channels were used in the propagation. Resonances converging to the first 4 excited Landau 
thresholds are shown in figure [7J The resonance structure is that of strongly perturbed Rydberg resonances converging 
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FIG. 7: Photoionization cross section of helium in arbitrary units versus the energy measured relative to the field free ionization 
threshold. Excitation is from the ground state in a field strength of 4700T using linearly polarized light. 

to the individual thresholds. 

The spectra in both cases, when calculated over an extended energy range, repeat the patterns shown in figures [6] 
and 



VI. CONCLUSIONS 



We have presented a detailed description of a method to evaluate the photoionization cross section of an atom in an 
external magnetic field. By combining R-matrix propagation with local adiabatic bases we have shown that is possible 
to calculate the cross section over an extended energy range for a range of field strengths. We have calculated cross 
sections for a range of atoms in both laboratory and astrophysical field strengths to illustrate the generality of the 
method. In addition, for given values of B, m and tt z , the spectra of all atoms of interest can be calculated without 
the need for any additional propagations. By using MQDT one is able to calculate quickly all of the resonances in 
the spectrum and to analyze some of their main characteristics. Partial cross sections to individual Landau levels 
are evaluated when calculating the total cross section enabling one to calculate their distributions. The method 
can be used to calculate the large amounts of data needed for such problems as stellar opacities or for calculating 
re-combination rates at low temperatures for an atom in a magnetic field. 



VII. APPENDIX 



In constructing the matrix representation of the Hamiltonian plus Bloch operator in each sector one has to construct 
matrices for the operators 

_ 1 d 2 + _ i 

2dr 2 + 2r 2 r W 

and r 2 with the radial basis set of shifted Legendre polynomials. 

For an arbitrary sector with inner radius r = a and outer radius r = b the coordinate r can be rescaled to the 
coordinate u such that 

o / / J, i „ \ \ 

(39) 
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The limits of the integrals in r over the limited range a — » b therefore become — 1 — ► 1 in the coordinate u. The 
orthogonal radial functions are thus 



Mr) = y yr^n-xfo)- (40) 

where P n is the Legendre polynomial of order n. The kinetic term can be evaluated by integration by parts and using 
the relation (Copson page 282) 

+1 dP n (u) dP m (u) 

— ; ; du = n(n + 1) it m — n even 

du du 

= if m-n odd (41) 

for m > n, hence 



1 f b dfndf?n d = 1 V(2Ti-l)(2 m-l) 
2J a dr dr r 2 (b - a) 



dr = — — — n(n — 1) if m — n even 



if m-n odd (42) 



where m > n. 

The potential term I nm is given by 



' q AW J /m(r)dr = ^ (-1) ^ du. (43) 



Using the relation 

Pn {p^jP-m (^) 



dx = 2P n (z)Q m (z) (44) 

1 z — x 

where m> n and Q m {z) is a Legendre function of the second kind [22] the potential term becomes 



✓ <*. - !)(*» - 1)^ (_2£4\ Q ^ (_^ay _ 



Since 



one gets finally 



b — a \ b — a J \ b — a 



P n (-z) = (-l) n P„(z) 
Q n (-z) = (-l) n+1 Q n (z) (46) 



= V(2n- 1 )(27n-l) _ A + a\ 

o — a \o — a J \b — a 



The centrifugal term requires the integrals J n 



1 - , V(2n-l)(2m-l) P„_ 1 (u)P m _ 1 (u) 



Jnm = / fn{r)^f m (r)dr = v v— - -/v-- — :/ / - ^f-t- du. (48) 

6-o 



U- 1 fcH 



The integral in this equation can be evaluated by differentiating eqn. (|44p to obtain 

^ 1 P " (x)Pm ^ } = _^("( n+ l)p„ +1 ( z )Q m ( z ) + ( m + l)P n ( z )Q m+1 ( z ) + 

(n + m + 2)zP„(z)Q m (z) J (49) 
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where the following relations have been used 



(z 2 -l)^jf± = (n+l)P n+1 (z)-(n+l)zP n (z) 
{z 2_ 1) dQ^z)_ = {m + 1)Qm+l{z) _( m + 1)zQm{z y (50) 

Substituting the appropriate value for z and using the relations in eqn. (I46j) the final result is 

Jmn . iB^S M ( . ir ^ iM + 

mP n _i(c)Q m (c) -c(n + rn)F n _ 1 (c)Q m _ 1 (c)^, (51) 

where 

c=^. (52) 
o — a 

To evaluate the integral involving r 2 , the recurrence relation 

(2n + l)uP n (u) = {n + l)P n+1 {u) + nP„_ x (u) (53) 

is used. This yields 

6 t ,s 2 , , s , (b-a) 2 ( n(n + l) 

fn{r)r 2 f m (r)dr = ± — — v ; = =5 mn+2 



4 V (2n + 1) V2n - l\/2n - 3 

2cn ( 2 n 2 (n — l) 2 

=o m „+i + c + 



v/^r^TV^TTT V (2n-l)(2n + l) (2n-l)(2n-3) 

2c(n-l) x (n-l)(n-2) \ 

-1 + 7^ , 7= ,„ = <W-2 ■ (54) 



V2n- lV2n-3 (2n - 3)V2n - V2ra - 5 

To evaluate equations eqn. (|47|) and eqn. (fBTj) the functions P n (c) and Q n (c) must be calculated. To calculate the 
Legendre polynomials P n {c) the standard recurrence relation 

Pn+l(c) = (2 "+ 1)C P n (c) - -jTPn-lW (55) 

n + 1 n + 1 

is used where the first two Legendre polynomials are given by Po(c) = 1 and Pi(c) = c. The method required to 
calculate Q n (c) depends on the value c. For c < 1 the same recurrence relation as eqn. (|55"T) can be used with 

Qo(c) = \ In ffzi) and Qi(c) = f In ffzij ~ 1- For the case c = > 1 a different method is used. For c > 1 the 

recurrence relation in eqn. (fS"5"|) should only be used for decreasing values of n. The first two values can be evaluated 
using the expression 

Q„(c) - ^Pn(c)ln (\^) - W n -i(c), (56) 



2 " w V 1 -c 

where 



2n — 1 „ , . 2n — 5 „ , . 2n — 9 „ 
^n-x(c) = ^P^( C ) + 5 ^P n _,( C ) + s ^P n _ B (c) + 



n -, 

V -P m . 1 {x)P n . m (x). (57) 



m 

m—l 



The first values could also be calculated using hyper-geometric functions via 



7T3 r(n + l) 1 „ A n 1 n 3 1, 

' (c) ^rtr|y^( 1 + 2'2 + 2 ; " + 2v,) 
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for |c|> 1, however in the calculations described in this paper the expression in eqn. (|56|) was used. The Legendre 
functions of the second kind can therefore be evaluated using the recurrence relation 

n + 2 In + 3 

Qn(c) = —-Qn+2(c) — cQ„ +1 (c). (59) 
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